Ab initio derivation of Hubbard models for cold atoms in optical lattices 
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We derive ab initio local Hubbard models for several optical lattice potentials of current interest, including 
the honeycomb and Kagome lattices, verifying their accuracy on each occasion by comparing the interpolated 
band structures against the originals. To achieve this, we calculate the maximally-localized generalized Wannier 
basis by implementing the steepest-descent algorithm of Marzari and Vanderbilt [N. Marzari and D. Vanderbilt, 
Phys. Rev. B 56, 12847 (1997)] directly in one and two dimensions. To avoid local minima we develop an 
initialization procedure that is both robust and requires no prior knowledge of the optimal Wannier basis. The 
MATLAB code that implements our full procedure is freely available online at http : / /ccpf orge . cse . rl . 
ac .uk/gf /project/mlgws/ 
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Atoms loaded into periodic optical potentials (T] |2) can 
be sufficiently cold to only occupy a small number of low- 
est energy Bloch bands. The interaction between two atoms 
occupying the same potential well can be large so that they 
form a paradigm test-bed for studying the physics of strongly 
correlated quantum lattice models. To derive accurate mi- 
croscopic models it is desirable to express the state of the 
atoms in terms of a basis of highly-localized single-particle 
states, given by some unitary transformation of the Bloch 
states forming the lowest energy bands. The reasons for this 
are two-fold. First, the occupations of localized basis states 
are measurable through high-resolution imaging |3). Second, 
the Hamiltonian rewritten in terms of localized basis states is 
typically a Hubbard model dominated by a few local terms. 
Together these two points justify the simulation of local Hub- 
bard models, used to describe many phenomena in condensed 
matter, using cold atoms in optical lattices [1|. In this arti- 
cle, we develop a procedure to systematically find a set of 
highly-localized basis states and thereby derive ab initio the 
parameters of a Hubbard model realized using cold atoms and 
an optical lattice. 

Only in simple cases, e.g., a lattice potential that is or- 
thogonal JU or leads to an isolated lowest Bloch band Bl . 
have the parameters of Hubbard models realized by cold 
atoms in optical lattices been derived using a basis of lo- 
calized single-particle states. The single-particle states used 
are Fourier transforms of the Bloch states, called Wannier 
states [6). For more complicated optical-lattice potentials, 
Hubbard parameters have been estimated rather than derived 
from first principles: on-site interaction Hubbard parameters 
have been estimated by using Gaussians centered at lattice 
minima as approximations to the single-particle states, and 
nearest-neighbor hopping parameters found by fitting a tight- 
binding form to the energy structure of the bands, without 
a rigorous justification of the tight-binding assumption (see 
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e.g. Refs. J2] [7] O). The approach we take here improves 
upon such calculations in two ways. We use a class of single- 
particle states that generalize the Wannier states and can thus 
be more localized. Also, our procedure calculates Hubbard 
parameters from first principles, without approximation, and 
provides a quantitative justification of neglected terms. The 
necessity of such improvements has recently been noted in 
the literature (9). 

Our procedure is an adaptation of several others already in 
use in solid-state physics. Specifically, we take as our starting 
point an algorithm developed by Marzari and Vanderbilt iflOl . 
They consider a basis of generalized Wannier states; Fourier 
transforms of inter-band mixtures of Bloch states. Choosing 
some initial basis, a steepest-descent minimization algorithm 
is used to iteratively generate another set of generalized Wan- 
nier states with a smaller spatial spread. The desired end-point 
of these iterations is the basis corresponding to the global min- 
imum of the spread, the so-called maximally-localized gen- 
eralized Wannier states (see Ref. ifTTl and references within 
for a review on the topic). Once this optimal basis is found, 
the parameters of the corresponding Hubbard model are easily 
calculated. 

The currently available software packages ifTZll that im- 
plement the steepest-descent minimization algorithm operate 
in three dimensions. For use with optical lattice potentials, 
which are often effectively one or two-dimensional, we have 
implemented the algorithm directly in these lower dimen- 
sional spaces, as well as in three dimensions. We find that 
for the optical-lattice potentials considered here, our imple- 
mentation in conjunction with commonly used initialization 
procedures (e.g. that described in ifTUl ) typically fails to con- 
verge to the global minimum of the spread and instead be- 
comes trapped in a local minimum; the maximally-localized 
generalized Wannier states are not obtained. 

Therefore, our algorithmic contribution is a new initializa- 
tion procedure for the Marzari and Vanderbilt steepest-descent 
algorithm. Our initialization procedure has an additional ben- 
efit in that it requires no knowledge of the optimal Wannier 
states, e.g., their location or approximate form, and therefore 
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requires no input beyond specifying the lattice potential. The 
initialization procedure is split into two parts, each minimiz- 
ing the inter- and intra-band contributions to the spread of 
the generalized Wannier states, respectively. The former is 
a method for minimizing the spread in the case of a single 
band ifTUl . The latter relates to a procedure devised by Souza, 
Mazari and Vanderbilt to optimally disentangle a subset of 
bands from a group of degenerate bands [ 13 1. Our whole pro- 
cedure, taking the lattice potential as input, and outputting the 
maximally localized Wannier states and Hubbard parameters, 
is combined into a single MATLAB routine. We have made 
this code freely available online fl4l . 

Note that while in the last stages of preparing this article we 
became aware of a very recent article 1 15 ] in which the authors 
use a different procedure to compute the maximally-localized 
generalized Wannier states and justify a local Hubbard model 
for bosons in the two-dimensional honeycomb potential. 

The remainder of the article is organized as follows. In 
Sec. Uwe discuss the derivation of Hubbard models for cold 
atoms in optical lattices, introducing generalized Wannier 
states as a basis for this derivation and outlining the problem 
of finding the states with minimum combined spread. Our 
approach for obtaining the maximally-localized basis is then 
described in Sec. [H] We include an outline of Marzari and 
Vanderbilt's steepest-descent algorithm, discuss the steps of 
our initialization procedure and then summarize how we com- 
bine these elements. In Sec.[ni]we derive Hubbard models for 
bosons in several optical-lattice potentials, first in one dimen- 
sion then in two, verifying the accuracy of our calculations on 
each occasion. Finally, we conclude in Sec. IV before pre- 
senting computational details in the appendices. 
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obeying the translational equivalence 

i4(r) = w R ,(r+R'-R). 
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Here R is a direct lattice vector for which V(r + R) = V(r) is 
satisfied, and which indicates the lattice site where w R (r) is 
localized, relative to some origin. The integer n is commonly 
called the band number, although as we shall see shortly it 
will index modes which may comprise of mixtures of several 
bands. An atom occupying the mode w' R (r) is often said to be 
in the n-th excited state or mode of lattice site R. 
The expansion thus takes the form 
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where Z? R annihilates a boson in mode VfS(r), such that the 
Hamiltonian H may be re-expressed as 

mn RR' 

1 1 v v jT mn °p h m ^h n ^h° y 

+ 7 U RR'R"R"' £, R "R'W" 

Z m™pRR'R"R"' 

with hopping and interaction parameters 



I. OBJECTIVE 
A. Hubbard models for atoms in optical lattices 

To begin, we outline the typical approach to deriving Hub- 
bard models for ultracold atoms with mass /j in an opti- 
cal lattice. For simplicity we assume the atoms to be spin- 
less bosons; extensions to fermionic atoms, multi-component 
gases including Bose-Fermi mixtures with different lattice po- 
tentials, atom-molecular interactions and finite-range interac- 
tions are straightforward JT]. 

Standing waves of laser light, tuned out of resonance, exert 
a spatially-periodic AC Stark shift on the ground internal state 
of the bosons. For sufficiently low atom energies and densities 
p the interactions between the atoms are well-approximated 
by a contact interaction of strength g. The effective Hamilto- 
nian is then of the form |4] 

H = j dr4' t (r)M'(r) + |y' dr ^(r)^(r)4f(r)4f(r). 

Here W annihilates a boson of mass /j and the single-particle 
Hamiltonian is h = — /i 2 V 2 /2 J u + V(r), where V(r) is the lat- 
tice potential induced by the AC Stark shift. 

We expand the field operators in terms of a complete ba- 
sis of orthonormal mode functions w R (r), corresponding to 



U m>W» = sf dr<(r)^(r)^(r)<„,(r). 

Due to Eq. (JTJ, these parameters are invariant under a simul- 
taneous translation in the direct lattice vectors that label them. 

The Hamiltonian simplifies in two ways. First, for suffi- 
ciently small kinetic E^ n and interaction energies Ei nt « pg, 
we can ignore all but some number / of the bands. Second, 
Vfn(r) are chosen such that they are well localized, meaning 
that the tj™, and ^rd/rwrw corresponding to hopping or inter- 
action between distant states are negligible. This leaves the 
Hubbard model 

J _ 

HhM = " £ £ f RR'^R b K> 

mn=\ (RR') 

1 J 

+- v v jf'"""p h mf h'^h° y 

+ 7 L i-i u rr'r"r'""r "r'W" 

z mnop=l (RR'R"R'"> 

where the angular brackets indicate that the sum is restricted 
to local terms, e.g., same-site, nearest-neighbor, or next- 
nearest-neighbor etc. The range of the terms that need to be 
kept will depend on how local the w R (r) can be, which in turn 
is dependent on the form of the potential V(r). 
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B. Generalized Wannier states 

We now turn our attention to the choice of wavefunctions 
Wg(r) used in the above procedure. A complete basis of or- 

thonormal functions is provided by the Bloch states |Vm )> 
corresponding to eigenfunctions of h 

¥ i k) (r)=e lk - r «i k) (r), 

where u)„ (r) are cell-periodic functions 11161 . The Bloch 
states of a given band m are uniquely labeled by a wave-vector 
k that runs over the first Brillouin zone of the reciprocal lat- 
tice. Any band m with energies 

E^ = {^\h\^}, 

(k) 

satisfying E„ t T^in, Ei nt for all k will not contribute to the 
physics and may be ignored. For all optical lattice potentials 
we consider here it is possible to focus solely on a small num- 
ber 7 of the lowest-energy bands which may be degenerate 
amongst themselves but are separated in energy from the oth- 
ers. 

To describe local interactions within this /-band subspace, 
a good choice of basis are states of the form 

\Rn) = 7^o[ dke-' k ' R £ f/^| ¥ i k) ), (2) 
{2%)" Jbz ^ 

where T is the volume of the primitive cell of the D- 
dimensional direct lattice, and C/W [ s a unitary matrix that 
mixes the Bloch bands. In the case that U^ k ' is diagonal, 
i.e., there is no band mixing, these states are exactly those 
first considered by Wannier J6). Thus the states appearing 
in Eq. Q are commonly referred to as generalized Wannier 
states. 

The separation in energy of the 7 lowest bands from the oth- 
ers ensures that some states |R«) exist with mode functions 
Wg(r) that are exponentially localized at lattice site R in co- 
ordinate space 02H2T). This exponential localization occurs 
if and only if the Bloch superpositions 

£^ k) | ¥ , (k) ), (3) 

m— 1 

are analytic (infinitely differentiable) in k across the whole 
Brillouin zone 11221 . This is a rigorous way of saying that only 
smoothed-out Bloch superpositions will lead to localization 
when Fourier transformed. When there are no degeneracies 
between bands can one simply use the phases of elements 
of a diagonal t/' k ' (representing the freedom in the phase 

of each |\|/m^)) to ensure the smoothness of the Bloch states 

|\j/i k ^). Hence simple Wannier states provide an exponen- 
tially localized basis in such cases. However, this is no longer 
the case when degeneracies and crossings in the band struc- 
ture lead to non-analytic |v|/j k '). In this situation band mix- 
ing and therefore a non-diagonal are required to obtain 
smooth Bloch superpositions and an exponentially-localized 



basis. The 'only if case highlights the importance of the gen- 
eralization of Wannier states to include non-diagonal t/W. 
Even when exponential localization is possible using simple 
Wannier states, generalized Wannier states may still signifi- 
cantly improve the localization. We will give examples of this 
inSecHn] 

C. Maximally-localized generalized Wannier states 

Generalized Wannier states therefore have the potential to 
provide a well-localized basis for the derivation of a Hub- 
bard model. However, generalized Wannier states are highly 
non-unique and so it remains to find and choose a single 
exponentially-localized basis. 

Several criteria have been proposed as a means of select- 
ing a specific basis of generalized Wannier states fT0ll2"3"ll24l . 
Here, following Ref. |10|, we seek the generalized Wannier 
states with a minimal combined spatial variance, henceforth 
called spread, defined as 

g= EL [(o«|f 2 |o«)-(o«|f|o«) 2 ] 

= LLi[(* 2 )n-rl] =L J n=i n n . (4) 

The minimizing states are called the maximally-localized gen- 
eralized Wannier states. 

It is known that the maximally-localized generalized Wan- 
nier states are indeed exponentially localized [ 17-19, 25ll26l . 
They do not necessarily provide the optimal approximation to 
H when restricting the number of terms kept in Z/hm but can 
be expected to be close to optimal. Finding a set of general- 
ized Wannier states which makes the Hubbard model approx- 
imation optimal is challenging. Hence minimizing the spread 
of the generalized Wannier functions is both an effective and 
practical choice with the added benefit of having a straight- 
forward interpretation in terms of a particle occupying a spe- 
cific lattice site. Note that it has been hypothesized that the 
maximally-localized generalized Wannier states always cor- 
respond to real functions, up to a global phase, when dealing 
with an isolated group of 7 bands EOl l26l [271 . All of our 
calculations support this hypothesis. 

II. METHOD 

To obtain maximally localized Wannier states we must first 

calculate the band structure and Bloch states |v|/,^ ). Such 
calculations are well-understood and we include our proce- 
dure here only for completeness and to introduce notation re- 
quired later. Then we must find the gauge such that the 
resulting generalized Wannier states |Rn), defined by Eq. 0, 
minimize the spread, defined by Eq. Q. This usually consists 
of several steps: initially, we choose the gauge = lj to 
be the 7x7 identity matrix. Then the gauge is transformed it- 
eratively t/< k )y( k ) according to a unitary V^ k '. These 

transformations accumulate until they converge to the desired 
gauge U^ k \ corresponding to the minimum spread. From this, 
the maximally-localized generalized Wannier states may be 
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calculated together with the hopping and interaction parame- 
ters. 

In the strategy devised by Mazari and Vanderbilt 1 10 1 an ini- 
tial unitary V' k ) is constructed via projections of J localized 
trial orbitals onto the Bloch states. This transformation leads 
to a gauge U 1 ^ = corresponding to an analytic set of 
Bloch superpositions |tjirj, ) and thus exponentially-localized 
generalized Wannier states |Rn). Subsequently, other uni- 
tary transformations V( k > are iteratively applied as part of a 
steepest-descent algorithm, in the hope that the cumulative 
gauge converges towards the spread-minimizing gauge. 

Unfortunately when applying this strategy to investigate 
common optical-lattice potentials, we found that the gauge 
corresponding to the maximally-localized generalized Wan- 
nier states is rarely obtained. Instead the spread often con- 
verges to some non-global minimum. Therefore we adopt 
a different initialization procedure, to precede the same 
steepest-descent algorithm. Contrastingly, we find that, for 
the optical-lattice potentials considered, our strategy con- 
verges quite consistently to the global minimum, and thus the 
maximally-localized generalized Wannier states are reliably 
obtained. While proving convergence to the global minimum 
is difficult, we test numerically both that the same solution is 
found when starting from two Bloch states differing by ran- 
dom permutations of the bands at each k, and that the final 
generalized Wannier states are real. 

In this section we begin by outlining the band structure 
calculation, and the representation of quantities, such as the 
spread, in reciprocal space. We then briefly describe Mazari 
and Vanderbilt's steepest-descent algorithm, before outlining 
two methods we will use as our initialization procedure. To 
end the section, we describe our full procedure for calculating 
the maximally-localized generalized Wannier states. 



A. Band structure 

We work with the Fourier space representation of the po- 
tential and cell-periodic functions 



(k,G) 
Cm 



7= f drV(r)e~ 
/Tipc 
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FIG. 1: (Color online) Representing the problem in k-space. (a) The 
reciprocal lattice points for a two-dimensional oblique lattice. The 
black circle has a radius of G ma x and is centered on T. Fourier com- 
ponents corresponding to reciprocal lattice points within the circle 
(red points) are included in the truncated basis set, while those out- 
side (blue points) are not. (b) A mesh of wave-vectors k for a hexag- 
onal two-dimensional lattice used to interpolate values of functions 
of k over the Brillouin zone, which is shown by the black hexagon. 



magnitudes |G| less than G max , as shown in Fig.[TJa). Then we 
only solve Eq. |5]) for a D-dimensional uniform discrete mesh 
of Vf wave-vectors k = G/M contained within some primi- 
tive cell of the reciprocal lattice, and interpolate between these 
wave-vectors. As shown in Fig. [TJb), this primitive cell need 
not be the first Brillouin zone since the corresponding Bloch 
states are invariant when translated by a reciprocal lattice vec- 
tor into the Brillouin zone. For each k in the mesh, solving the 
set of Eqs. |5]) then reduces to an eigenvalue problem. The jus- 
tification of the truncation and mesh discretization, as well as 
the values of G max and M we use are discussed in AppendixfA] 
Note that various symmetries guarantee certain properties 
of the coefficients Cm ' G ' [28 1. For a given m and k, time- 
reversal symmetry implies that we may choose c m k ' G ' — 
Cm' G '*. The addition of inversion symmetry allows us to set 
Cm' G ' and v' G ' as real up to a common phase factor and en- 
sures y( G ' = y(~ G ). This implies that, when both symmetries 
are present, we may both reduce our mesh of reciprocal lattice 
vectors by nearly half, as the coefficients for G may be in- 
ferred from those for G, and restrict all quantities in the eigen- 
value equation <|3j to be real, thus speeding up computations 
for each k. 



where the integral is over a primitive cell of the direct lattice. 

In this representation the single-particle Schrodinger equation 

fi 
h\vf), 



= Em^^^m 1 ) may be written as 



2M iff £7 

The full band structure is obtained by solving this equation 
for all wave- vectors k in the Brillouin zone and all reciprocal 
lattice vectors G fl6l . 

To make this calculation tractable on a computer, we firstly 
truncate the Fourier expansions to include some finite number 
N of terms, corresponding to reciprocal lattice vectors G with 



B. Contributions to the spread 

Following Mazari and Vanderbilt, the spread £2 = Q.j + £2 
can be conveniently decomposed into two positive definite 
parts, £2/ and £2. The latter depends on the choice of gauge 
t/M appearing in Eq. ^2), while the former does not. The 
gauge-independent part £2/ depends only on the smoothness in 
k-space of the underlying manifold of Bloch states, while the 
gauge-dependent part £2 depends on the additional smooth- 
ing achieved by applying phases to and mixing the Bloch 
states. In preparation for what follows, it is useful to further 
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decompose Q.j = ii/.o + Qi,od and Q. = Q.o + CIod into band- 
diagonal and band-off-diagonal terms. For the diagonal terms 
Q.j = £i" D and Q.£> = Q,' D it makes sense to break them 
down into positive-definite contributions from each band. The 
decomposition is expressed neatly as 



It is clear then that all the information about the spread is 
contained in the matrix elements 

(k,b) _ (k)* v (k,G)* (k+b,G) (k) 
w*mn — 2^ P m 2^ P n ' 

op G 

which are the truncated Fourier representation of the overlap 
i k+b) ), where similarly to Eq. §, 



Our minimization method will of course leave £2/ invariant, 
while minimizing CI. The Mazari and Vanderbilt steepest- 
descent algorithm iteratively minimizes Cl directly, while our 
initialization procedure is divided into two stages, one which 
reduces Hod and another which minimizes dp- 

In terms of generalized Wannier states, the contributions to 
the spread are written 

a n ip = (0«|f 2 |0«)-£|(0«|f|R«)| 2 , (6a) 

R 

£1,,od = -£££|(0m|f|R«)| 2 , (6b) 

n m^=n R 

a D = £ |(o»|f|R»)| 2 , (6c) 

R^O 

^od = £££|(0m|f|R«)| 2 . (6d) 

n m^n R 

Again, for computational tractability, we move to the trun- 
cated Fourier representation with a discretized mesh. In this, 
all integrals over the Brillouin zone are replaced by summa- 
tions over the mesh, 



{2%) D 



dk 



BZ 
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and gradients represented by finite differences (the gradients 
in reciprocal space arise from moments in position space). We 
use the finite-difference expressions recommended by Marzari 
and Vanderbilt [ 10 1, which have the property of transforming 
correctly under translations of the generalized Wannier states 
by a direct lattice vector. In this way, contributions to the 
spread are re-expressed as 



i«n = ii/&v>, 



(8) 



with \um ) the state associated with periodic function M^'(r). 

These elements are initialized to Mmn = ^GcJ' G '*cf i+l, ' G ' 
when t/( k ' = 1/. Then under a gauge transformation U^> — > 
j/( k )y( k ) they undergo the computationally simple transfor- 
mation -> y(k)t M (k.b)y(k+b)^ 



C. Minimizing total spread 



The gradient r' k ' = dD./dW^ k \ embodying the change in 



with 



spread due to a gauge transformation 
dW( k ) an infinitesimal anti-Hermitian matrix, can be effi- 
ciently calculated from the matrices M' k b ' (see Ref. iflOl for 
details). The steepest-descent approach, as used in Refs. ifTUl 

[T3l , then implements the gauge transformation V^ k ' = c dw k) 
with dW^ = — eT^ and e a small positive number. These 
steps are repeated until convergence is achieved. 

To a large extent, the steepest-descent algorithm is only 
as good as its initialization procedure, since starting from an 
arbitrary set of generalized Wannier states, the algorithm is 
likely to drive the set towards one of the many local minima in 
the spread, rather than the global minimum. We do not discuss 
here the commonly used projection-based initialization proce- 
dure (see Ref. ifTTI for information on this) that we found to 
struggle for optical-lattice potentials. Instead we now discuss 
two other approaches for reducing the spread, which together 
will form the initialization procedure we use successfully for 
optical-lattice potentials. 
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(7a) 
(7b) 



^CO b (lm[lnMl k ' b) ]+b-r„) 2 , (7c) 

1 k,b V ' 



(7d) 



Here the vectors b connect each wave-vector k to its nearest- 
neighbors, C0b are factors that depend on the geometry of the 
mesh [29 1, and 



^£co b bIm[lnM, 



k.b 



D. Reducing inter-band spread 

We break down the task of finding the maximally-localized 
Wannier states into two stages. The first stage is to mix 
the bands to create a new set of pseudo-bands from which 
a maximally-localized ordinary Wannier states calculation is 
optimal (i.e. leading to the smallest possible spread). The 
second stage is to calculate the maximally-localized ordinary 
Wannier states using these pre-mixed bands as a starting point. 
The first stage corresponds to minimizing the off-diagonal 
term £2od, and the second to minimizing the diagonal term 
Q.D- Our initialization procedure is split accordingly: first we 
reduce (but not necessarily minimize) Q.od, as described in 
this subsection; second we minimize £2q, as described in the 
next subsection. 
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Our first goal is then to reduce the band-off-diagonal term 
CIod, which is equivalent to reducing £2/£>. This equiva- 
lence is clear from the interpretation above, that reducing 
Q.OD corresponds to optimizing the bands from which to per- 
form a maximally-localized ordinary Wannier states calcula- 
tion. Mathematically, it follows from observing that the band- 
off-diagonal parts of the gauge-invariant spread £2/ and the 
gauge-dependent spread £2 are the negative of each other (see 
Eqs. (J6]?), <[6ji) and (jT]?), |7}i)): reducing CIod is achieved 
by increasing the band-off-diagonal part £2/.od of the gauge- 
invariant spread or, equivalently, reducing its diagonal part 

&I,D- 

To reduce £2/ £>, we use a method devised by Souza, Mazari 
and Vanderbilt [13]. For K degenerate bands, their minimizes 
the contributions to £2/ from a subset K' < K bands obtained 
through a unitary mixing of these bands. The aim of this ap- 
proach is then to construct the K' bands with the smoothest 
k-space such that they provide the optimal set of K' bands 
from which to construct localized generalized Wannier states 
(optimal in the sense of having the smallest possible gauge- 
invariant contribution to the spread). 

We use their approach to reduce £2/o in the following way: 
First, we use the Souza et al. method to minimize £2) D and 
therefore construct, from the K = J bands, a single (K' = 1) 
band whose smoothness in k-space is optimum for construct- 
ing a localized Wannier state. Then, keeping this band fixed, 
we use the Souza et al. method again to minimize £lj D and 
construct from the K = J— 1 remaining bands a single (K' = 1) 
band that is optimum for constructing a localized Wannier 
state. This is repeated in a similar fashion to obtain a third, 
fourth, etc. band until finally we use the Souza et al. method 
to minimize £2^ z) 1 and construct an optimized (/— l)-th band 
out of the two remaining bands, with all lower bands fixed. 
Our approach therefore consists of J — 1 applications of the 
Souza et al. method, in each case optimally extracting a sin- 
gle (K' = 1) band from K = J,J-1,...,2 others. 

Note that this does not necessarily minimize £2/ j> and there- 
fore Hoc but we find that following this procedure Q.od is 
very small. Details of the Souza et al. method and our use 
of it can be found in Ref. [ 1 3 1 and Appendix [B] respectively. 
Here we simply note that the Souza et al. method proceeds via 
several iterations, each of which applies a transformation 
over all k-space that would have minimized the contribution 
to the spread from any given point in k-space had there been 
no transformation applied at the other points in k-space. The 
desired gauge must be left unchanged by such an iteration and 
thus it is a possible point of convergence. To protect against 
false convergences, we initialize the whole procedure above 
by applying a transformation V^ k ', where at each k we take a 
J xJ identity matrix and randomly permute its rows. We find 
that in practice, following this initialization, the desired gauge 
is nearly always obtained. 



E. Reducing intra-band spread 

We next present a method that reduces the intra-band con- 
tribution £2o to the spread, while leaving the inter-band con- 



tribution £2<3£> invariant. To ensure this invariance, in this sec- 
tion we restrict ourselves to gauge transformations V* k ) that 
are diagonal, i.e., while we allow changes to the phases of the 

Bloch superpositions |vj/i k '), we do not allow any transforma- 
tions that further mix the bands. Hence the task splits into 
J independent parts, each to reduce £2£, by applying phases 

(k) 

to the Bloch superpositions |\j/„ ). One may interpret this as 
constructing the maximally-localized ordinary Wannier states 

from a set of bands comprising the mixed Bloch states l 1 ^'). 
Such single-band tasks are usually described in terms of 

the Berry connection A< k > = i(u ( n ] |V( k ) \u { n ] ) GUI. Integrals 
(Berry phases $c) of the Berry connection around closed 
paths C in the Brillouin zone are invariant under changes to 

the phases of the Bloch states |v/i k '). This implies that f„, 
equal to the average value of A^ k ^ across the Brillouin zone, 
is also invariant [31 1. A further invariant quantity is given 
by B = V x A' k \ called the Berry curvature. Local values 
of A' k ', however, depend on the phases of the Bloch states, 
which determine the phase-dependent part of the spread £2^,. 
It is known that this spread is minimized when the divergence 
of the connection vanishes, V • A^ k ' = 0, and the minimum 
possible spread depends only on the Berry curvature B. 

In particular, if B = 0, then the minimum possible spread 
£2£, is zero. It follows that all Berry phases are zero and it is 
possible to smooth A' k ) such that it is uniform, at which point 
£2q = 0. To smooth the connection A' k \ we use a progres- 
sive phase update method: it consists of taking a succession 
of closed loops through the Brillouin zone, and, for each, al- 
tering the Bloch phases at points along the loop such that the 
projections of A^ along it are constant. This constant value is 
fixed by their integral around the loop, the Berry phase, which 
is invariant. Adjusting the phases in this way for several loops, 
given in Appendix [C] will result in a flattened connection, if 
possible. 

For non-zero B, absolute uniformity of the connection A( k ) 
is not possible. However, in an attempt to suppress the di- 
vergence of the connection and therefore approach the min- 
imum spread £2^,, we still choose to smooth out A' k ' using 
the progressive phase update method and find this greatly re- 
duces £2£,. To achieve the minimum, we follow the progres- 
sive phase updates with the steepest-descent minimization al- 



gorithm of Marzari and Vanderbilt (cf. Sec. EC I, when only 
terms corresponding to £l" D contribute to the gradient r' k ' . 

A particular case of interest is a system with inversion sym- 
metry and a current gauge that is diagonal, i.e., the bands 
have not been mixed. As a result of the symmetry, the mode 

functions u^ 1 (r) °< w„ (r) must be real up to a global phase, 
at which point the Berry curvature and thus the spread £2£, 
vanishes. Note that optical-lattice potentials usually possess 
inversion symmetry since this is inherited from the lasers that 
created them; superlattice techniques are required to break 
this. 

We found that even if f/' k ' is not diagonal, e.g., after the 
inter-band spread is reduced, the output of the disentangling 
procedure often still had zero Berry curvature for each band 
and the progressive phase update method reduced Q" D to zero. 
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1 . Calculate band structure 



2. (a) Apply progressive phase update to each band 
(b) Apply steepest-descent minimization to each band 



Ordinary Wannier states 



Generalized Wannier states 




3. Apply disentangling procedure to composite group 
I 

4. (a) Apply progressive phase update to each band 
(b) Apply steepest-descent minimization to each band 

I 

5. Apply full steepest-descent 



1_ 



6. Calculate the maximally localized Wannier states on a real space grid, 
and compute the Hubbard parameters 
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FIG. 3: (Color online) One-dimensional superlattice. (a) The con- 
figuration of lasers red-detuned from wavelength A, to produce the 
superlattice potential, (b) The potential over the unit cell, for s = 
(blue solid line), s = 0.5 (red dotted line) and 5=1 (green dashed 
line), (c) The band-structure corresponding to the potentials in (b). 



HI. RESULTS 



FIG. 2: Flow diagram of how our software package calculates the 
maximally-localized generalized Wannier states. First, it calculates 
the band structure, and from this computes the matrix elements 
c m ' c n + ' • Second, it minimizes £1'^ for each band 
as far as possible without mixing the bands. Third, it reduces Hod- 
Fourth, it again minimizes Q.'^ for each band as far as possible with- 
out further mixing the bands. Fifth, it minimizes the total £1 to its 
global minimum via steepest-descent minimization. Last, we com- 
pute the Hubbard parameters and construct the maximally localized 
generalized Wannier states. 



Specifically, this occurred whenever the degeneracies in our J- 
band subspace were a result of purely geometric symmetries. 
We hypothesize this is a general feature, also hinted at in the 
results ofRefs. |20l |26l |271 . 



F. Full procedure 

Having described the elements of our computational ap- 
proach, we now describe how they are pieced together. The 
full procedure for calculating the maximally-localized gener- 
alized Wannier states is shown in Fig. [2] 

First we calculate the band structure. Then we minimize the 
intra-band spread via a progressive phase update, followed by 
a restricted version of the steepest-descent method (for po- 
tentials with inversion symmetry, the steepest-descent part is 
unnecessary). At this point we are at the gauge corresponding 
to the maximally-localized ordinary Wannier states. We then 
reduce the inter-band spread using the method adapted from 
Souza et al. [ 13 1. The inter-band spread reduction usually has 
the side-effect of increasing the intra-band spread slightly, so 
we again apply a progressive phase update, followed by the 
restricted version of the steepest-descent method to minimize 
the intra-band spread. The above forms our initialization pro- 
cedure. If this has not already found the global minimum of 
spread, we find that it is sufficiently close that the full steepest- 
descent algorithm iflOl returns the maximally-localized gener- 
alized Wannier states with a close to perfect success rate. 



We next use the above procedure to derive ab initio the 
Hubbard Hamiltonians realized by bosons in a variety of 
optical-lattice potentials. The reasons are four-fold. First, 
we test the accuracy of our procedure. Second, we com- 
pare the procedure against others, e.g., methods using ordi- 
nary rather than generalized Wannier states. Third, we demon- 
strate that local Hubbard Hamiltonians can be justified for sev- 
eral experimentally-important cold-atom optical-lattice sys- 
tems. Fourth, we provide the relevant model parameters ac- 
curately in terms of well known control parameters like the 
laser intensity. 

For the testing, we use the hopping parameters from our de- 
rived Hubbard Hamiltonian to calculate an interpolated band- 
structure according to the tight-binding model. The legitimacy 
of our approximations can then be considered by comparing 
the interpolated band-structure to the original. Note that this 
only allows us to determine the accuracy of the hopping pa- 
rameters, not the interaction terms. Since we are unable to 
directly verify the accuracy of discarding interaction param- 
eters, we only discard those of magnitude equal to, or less 
than, that of the discarded hopping parameters for some typi- 
cal range of interaction strengths g <g = Ef{k D , where 



is the recoil energy, and X is the 'averaged' wavelength of the 
laser beams creating the optical lattice ll40l . 

We now obtain the maximally-localized Wannier states 
and nearest-neighbor Hubbard models for atoms in several 
one- and two-dimensional optical-lattice potentials. We leave 
three-dimensional potentials for a future presentation. 



A. One-dimensional systems 

To begin, we find the maximally-localized generalized 
Wannier states and related Hubbard parameters for bosons in 
a one-dimensional superlattice potential, given by 

V(x)=V [(l-s)sin 2 (27tx/A.)+.?sin 2 (47tx/A.)] . 
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(a) (b) (c) 




FIG. 4: (Color online) Maximally-localized Wannier states for the 
one-dimensional superlattice with Vo = 20£j? and s = 0.999. (a) The 
dotted green and dotted-dashed blue lines are the m = 1 , 2 general- 
ized Wannier states, (b) The solid red and dashed light blue lines are 
the n = 1,2 maximally localized ordinary Wannier states, (c) The 
spreads of the maximally localized ordinary and generalized Wan- 
nier states as a function of s (line type is the same as in (a) and (b)). 




0.5 g 1 0.5 j 1 0.5 s 1 



FIG. 5: (Color online) Hopping parameters for the one-dimensional 
superlattice. The lines show the magnitudes \t" m \ = \t^\ with 2|R| = 
jX. The black, blue, red, and green lines correspond to j = 0, 1 , 2, 3, 
respectively. The solid (dashed) lines correspond to the maximally- 
localized generalized (ordinary) Wannier states. 



Such a potential can be produced using two independent pairs 
of laser beams, each red-detuned from wavelength A, and at 
an angle to each other, as shown in Fig. [3|a). Their total in- 
tensities determine the potential depth Vo and their relative in- 
tensities determine the superlattice parameter < s < 1. The 
lattice parameter is X/2. This system has been experimen- 
tally realized in Refs. l32l[33l . and was proposed in Ref. P4l 
as a method for initialising a quantum register on a time- 
scale that is an order of magnitude smaller than the conven- 
tional quantum-freezing of a superfiuid to a Mott insulator 
state 1411351. 

The potentials for three values of s are shown in Fig.|3jb), 
where Vq = 20£r. For s = the potential is sinusoidal with a 
minimum at the center of each primitive cell. For s ^ there 
are two minima in each primitive cell, which move either side 
of the center. As s — > 1 the potential approaches a sinusoid 
with lattice parameter A/4. 

The band structures for the same parameters are shown in 
Fig. [3jc). For all < s < 1 the two lowest lying bands are 
well separated from the higher bands, and also are not degen- 
erate amongst themselves. Therefore, ordinary Wannier states 
will provide an exponentially-localized basis. We use this su- 
perlattice potential then to demonstrate that using generalized 
Wannier states can further localize the Wannier states even 
when there are no inter-band degeneracies. 

We expect the benefits of generalized over ordinary Wan- 
nier states to be most notable in the s — >• 1 limit, where the 




FIG. 6: (Color online) Interaction parameters for the one- 
dimensional superlattice. The lines show the magnitudes \UJ" l \ = 
I^OORrI °f interactions between two particles in bands m and n at 
sites separated by 2|R| = jX. The black, blue, red, and green lines 
correspond to j — 0, 1 , 2, 3, respectively. The solid (dashed) lines cor- 
respond to the maximally-localized generalized (ordinary) Wannier 
states. 



(a) (b) 




FIG. 7: (Color online) Accuracy of the one-dimensional superlattice 
Hubbard models, (a) Lowest and first excited bands for the super- 
lattice potential with Vo = 20Er and s = 0.999. The blue and red 
solid (cyan and magenta dashed) lines show the interpolated lowest 
and first excited bands, respectively, for a tight-binding model us- 
ing the maximally-localized generalized (ordinary) Wannier states. 
For maximally-localized generalized Wannier states, the interpolated 
bands are indistinguishable from the exact bands on this scale, (b) 
The standard deviation a between the exact bands and the interpo- 
lated bands as a function of the superlattice parameter s. The solid 
blue (dashed red) line is for maximally-localized generalized (ordi- 
nary) Wannier states. 

bands are close together. The maximally-localized general- 
ized Wannier states and ordinary Wannier states for the case 
s = 0.999 and J = 2 are presented in Figs. |4|a) and (b), respec- 
tively. It is clear from inspection that the generalized Wannier 
states are more localized. This improved localization occurs 
for even small s but is very significant for moderate or large 
s > 0.5, as is shown in Fig. |4jc), which plots the respective 
spreads as a function of s. 

Another way to see the effects of improved localization is to 
look at the magnitudes of the hopping and interaction param- 
eters. These are shown in Fig. [5] (hopping) and Fig. [^(interac- 
tion), as a function of s, for both the ordinary and generalized 
Wannier states. For large s > 0.5, non-local Hubbard parame- 
ters are significantly reduced when using generalized Wannier 
states. This comes at the expense of allowing inter-band hop- 
ping. 

From these values it is clear that using either ordinary or 
generalized Wannier states, a tight-binding Hamiltonian 

#hm = E t ( -4 n bp n j -srhf (b" i+l +b" J _ 1 ) 
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+ \uZ n b'fb'fV)b" ] + £ [-Cbfh) 

m= 1 ,m^n 




can be derived and justified from first principles. Here, for 
clarity, we have replaced the label R by the label j = 2\R\/X. 
The model derived using generalized Wannier states is more 
accurate, as we can demonstrate by comparing interpolated 
bands to the original. In Fig. [TJa) this is shown for both 
maximally-localized ordinary and generalized Wannier states, 
and superlattice parameter s = 0.999. The generalized Wan- 
nier states almost exactly reproduce the band structure, while 
there are significant deviations for the Hamiltonian derived us- 
ing ordinary Wannier states. In Fig.[7jb) we show the standard 
deviation,i.e., the average root mean squared error of the ener- 
gies averaged over the bands, between the interpolated bands 
and exact bands as a function of s. This demonstrates that 
the difference in accuracy between using ordinary and gener- 
alized Wannier states is appreciable for s > 0.5. In fact, we 
should have expected this from the non-sinusoidal nature of 
the bands for large s. A tight-binding model built from ordi- 
nary Wannier states can only ever result in a sinusoidal band 
structure. Generalized Wannier states and inter-band hopping 
they describe have no such restriction. 

These results confirm that for s > 0.5, the accuracy of 
the local model found using generalized Wannier states be- 
comes significantly better than that using ordinary Wannier 
states. The reason for this difference is that the two gener- 
alized Wannier states can each break the reflection symme- 
try in the primitive cell to localize around a different mini- 
mum (see Fig. Fig.[4ja)). Meanwhile the maximally-localized 
ordinary Wannier states cannot break this symmetry and in- 
stead are symmetric and antisymmetric combinations of two 
functions localized at each of the minima (see Fig. Fig. |4jb)). 
Aside from leading to more accurate local Hubbard models, 
the use of generalized Wannier states are more relevant for 
cold atom experiments. In such experiments, it is the presence 
of a particle at a position in space rather than the symmetry 
of its wavefunction that is measured through high-resolution 
imaging [3]. Thus a Hubbard model corresponding to atoms 
in spatially-separated sites is preferable to atoms in symmet- 
ric/antisymmetric superpositions. We similarly expect gener- 
alized Wannier states to be important for other lattices that 
possess more than one potential minimum per primitive cell. 



B. Two-dimensional systems 

The use of generalized Wannier states is paramount in two 
dimensions, as degeneracies in the lowest bands are likely to 
occur as a result of crystallographic point-group symmetries. 
Hence the maximally-localized ordinary Wannier states could 
fail to provide an exponentially localized basis due to the re- 
sulting non-analyticity of the bands. Further, we will see cases 
where the maximally-localized generalized Wannier states are 
not centered around inversion points of the lattice, and do not 
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FIG. 8: Hexagonal lattice, (a) The beam configuration for generat- 
ing the optical lattice. The three beams are blue-detuned from the 
wavelength X. (b) The lattice potential, with the white line marking 
the boundary of the Wigner-Seitz unit cell, (c) The band-structure 
for lattice depth Vq = WEr. The energies are displayed along the 
path through the Brillouin zone shown in the inset, (d) Similarly for 
V Q = 30£ R . 



l 



-l 




FIG. 9: Maximally-localized generalized Wannier states for the 
hexagonal lattice. The two lowest bands are shown for lattice depth 
Vo = 10£r. We have labeled the potential minima with equal hop- 
ping and interaction parameters from the 'home' minimum by j = 
0,1,2,3,4. 



share the symmetry of the lattice. In these cases, approxi- 
mating the states using Gaussian functions would lead to a 
particularly inaccurate estimate of the Hubbard parameters. 

To showcase our procedure we now calculate accurately 
and from first principles the maximally-localized generalized 
Wannier states and Hubbard parameters for atoms in an opti- 
cal lattice, with either hexagonal or Kagome geometries. Both 
potentials have multiple minima per primitive cell and lead to 
a degenerate set of lowest bands, thus representing a signifi- 
cant challenge using any other method. Both of these struc- 
tures also play an important role in condensed-matter physics, 
seee.g. Refs. lEESi^i. 



1. Hexagonal lattice 

Three blue-detuned beams of approximately equal wave- 
length X, shown in Fig. [8ja), generate a hexagonal optical- 
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FIG. 10: (Color online) Hopping and interaction parameters for the 



hexagonal lattice, (a) The magnitudes \tj 



ftnii 
l 0R\ 



of the hopping pa- 



rameters, as a function of lattice depth Vb, where the centers of \0m) 
and |Rn) are the y'-th smallest distance from each other (cf. Fig.|9|. 
The black solid, blue dotted, red dashed, green dot-dashed, and ma- 
genta dot-long dashed lines are for j = 0,1,2, 3,4 respectively, (b) 
Similarly for the magnitudes \Uj\ = |t/oORRl °f me interaction pa- 
rameters, (c) The total standard deviation a between the exact lowest 
bands and the interpolated tight-binding bands as a function of lattice 
depth. 



lattice potential, written as 

V{x,y) = — [3 + 2cos ) 

f3nx\ ( V3ny 

+4cos ( v irJ cos ( v — 

and plotted in Fig. [8jb). The potential exhibits two minima 
per unit cell, positioned at cell vertices that form a hexagonal 
(honeycomb) structure. The consequence of there being two 
potential minima per unit cell is that the two lowest bands are 
degenerate at the K points of the Brillouin zone, as shown in 
Figs.gc) and (d) for lattice depths V = 10£r and V = 30£ R , 
respectively. 

The maximally-localized generalized Wannier states for the 
two lowest bands, using 7 = 2, are shown in Fig. |9]for a lat- 
tice depth Vo = 10£r. Both states possess three-fold rotational 
symmetry about their centers and are images of one another 
through a rotation of 60° about the center of the Wigner-Seitz 
unit cell. As for the one-dimensional superlattice, both gen- 
eralized Wannier states are localized around a potential min- 
imum, rather than at the Wyckoff positions (centers of in- 
version). As a result ilo 7^ for the pair although the total 
spread is minimized and the cell-periodic superposed states 
(cf. Eq. |8]l) are real. Since inversion symmetry is broken it is 
clear that a Gaussian function would not adequately describe 
these Wannier states even in the deep lattice limit. 

The magnitudes of the hopping and interaction parame- 
ters for the two lowest bands are shown in Fig. [10] Since 
the maximally-localized generalized Wannier states are re- 
lated through a symmetry operation the parameters within 
each band are identical. We therefore label the parameters not 
by site and band, but by j, the rank of the distance between 
potential minima, as shown in Fig. [9] Parameters for j = 0,2 
are intra-band, while those for j = 1,3,4 are inter-band. We 
calculate these parameters up to a lattice depth of Vb = 200£r. 
We observe that for large Vb the significant parameters are the 
on-site interaction parameter Uq and the hopping parameter t \ . 
The interaction parameter U\ corresponding to the interaction 



FIG. 1 1 : (Color online) The Kagome lattice, (a) The beam config- 
uration for generating the optical-lattice potential. The projections 
of the beam wave-vectors are shown in the x-y and x-z planes. The 
beams are both red- and blue-detuned from the same wavelength X. 
(b) The resulting lattice potential, with the white line marking the 
boundary of the Wigner-Seitz unit cell, (c) The band-structure for 
lattice depth Vb = 2£r. The energies are displayed along the path 
through the Brillouin zone shown in the inset. 
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FIG. 12: (Color online) Maximally-localized generalized Wannier 
states for the Kagome lattice. The three lowest bands are shown 
for lattice depth Vo = 10£r. We have labeled the potential minima 
with equal hopping parameters from the 'home' minimum by j = 
0,1,2,3,4. 



of Wannier states in neighboring potential minima is also rel- 
atively large but is at least an order of magnitude less than t\ 
for Vb ^ 10£r and typical interaction strengths g s» g. With 
these observations, the Hamiltonian for the hexagonal optical 
lattice is accurately represented by 

Hum = - £ hSJBj + 'EkVot'Wihh (10) 

(U> 

where the sums are taken over potential minima (we have now 
dropped the band index and instead labeled the minima by the 
indices i and j), each with three nearest-neighbors, which we 
denote by the angled brackets. 

We once again insert the hopping parameters included in 
the Hamiltonian into a tight-binding model to recreate the 
single-particle band-structure. The standard deviation be- 
tween the interpolated bands and the exact bands is shown in 
Fig. 



10 c) as a function of lattice depth. This again decreases 



exponentially with lattice depth indicating the high accuracy 
of the model at all but shallow depths. 



2. Kagome lattice 

The Kagome lattice has received a large degree of interest in 
recent years because it leads to a highly-frustrated many-body 
Hamiltonian |7, 8 38 39 1. This lattice may be created using 
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deep lattice, to the Hamiltonian given in Eq. ( 10 1, where in- 
stead there are four nearest neighbors. Once again we can 
reassure ourselves of the accuracy of the derived Hamiltonian 



by looking at the standard deviation, shown in Fig. 13 c), be- 
tween the interpolated bands and the exact bands. This de- 
creases exponentially with lattice depth and is significantly 
smaller than the exact band-width, indicating good accuracy. 



FIG. 13: (Color online) Hopping and density-density interaction pa- 
rameters for the Kagome lattice, (a) The magnitudes \tj\ = \t^\ of 
the hopping parameters, as a function of lattice depth Vy, where the 
centers of |0m) and |Rn) are the j-th smallest distance from each 
other (cf. Fig. 1 12}. The black solid, blue dotted, red dashed, green 
dot-dashed, and magenta dot-long dashed lines are for j = 0, 1 , 2, 3,4 
respectively, (b) Similarly for the magnitudes \Uj\ = |t^)0RR I °f tne 
interaction parameters, (c) The total standard deviation a between 
the exact lowest bands and the interpolated tight-binding bands as a 
function of lattice depth. 



six lasers of approximate wavelength X, three of which are 
red-detuned and three of which are blue-detuned. The setup is 



shown schematically in Fig. 1 1 a) and the resulting Kagome 
potential 
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is shown in Fig. [TTJb). We scale the potential such that the 
full lattice depth Vb is the difference between the maximum 
and minimum of the potential V(x,y). The primitive unit cell 
possesses three potential minima, and the lowest three bands, 
shown in Fig. [TTJc) for a lattice depth Vb = 2£p>, are degen- 
erate; two of the bands are degenerate at the K points and are 
reminiscent of the lowest bands of the hexagonal lattice, while 
the highest energy band is almost flat and is degenerate at the 
r point. 

The three maximally-localized generalized Wannier states 
for the three lowest bands, using J = 3, are plotted in Fig. [12] 
for a lattice depth Vb = 10£r, and each is once again located at 
a potential minimum. Since the potential minima are located 
at Wyckoff positions, the generalized Wannier states possess 
inversion symmetry and £2d = 0. Each state is only two- 
fold symmetric under rotation in accordance with the point- 
symmetry of the Wyckoff position it is centered on. 

The states are images of each other through a rotation of 
120° about the center of a trimer. Because of this the mag- 



nitudes of the hopping parameters, plotted in Fig. 13 a), be- 
tween equivalent neighboring potential minima are equal, as 
was observed with the hexagonal lattice. Also similar to the 
hexagonal lattice, the parameters corresponding to hopping 
between adjacent minima decay almost exponentially, while 



the on-site interaction parameter dominates (see Fig. 13 b)). 
The nearest-neighbor interaction parameters are at least an or- 
der of magnitude smaller except at very low lattice depths. 
Once again, due to symmetry this leads us, for a sufficiently 



IV. CONCLUSIONS 

We have calculated, from first principles, the parameters 
of nearest-neighbour Hubbard models for several optical lat- 
tice potentials, including the honeycomb and Kagome po- 
tentials, demonstrating quantitatively for which lattice depths 
such models are accurate. Strongly-correlated phenomena 
probed in optical lattice experiments and quantum simulations 
depend delicately on the ratios of kinetic and interaction ener- 
gies. Therefore precisely determining them ab initio, as done 
here, is essential for diagnosing and interpreting such experi- 
mental results and for using optical lattices as quantum simu- 
lators. 

To perform our calculations we have developed a freely 
available software package lfT4l that, given an optical lat- 
tice potential, will efficiently calculate the corresponding 
maximally-localized generalized Wannier states without any 
prior-knowledge of their form in any spatial dimension. This 
will allow cold-atom researchers to easily and accurately de- 
termine Hubbard models realized by any laser setup. We hope 
that this tool will be useful for the optical-lattice community. 
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Appendix A: Accuracy of the truncated Fourier representation 
and discretized mesh 



The cut-off wave-vector corresponds to a maximum kinetic 
energy for the plane-wave components given by £ C ut-off — 
G^ ax /2M, and introduces a minimum spatial resolution A, m i n = 
27t/G max for describing real space functions in the system, 
namely, the potential, the Bloch states and the Wannier states. 
The minimum spatial resolution must be smaller than the spa- 
tial variations in these states in order for them to be accurately 
recreated using the truncated set of coefficients, therefore the 
cut-off energy must be at least as large as the highest energy 
Fourier component of the potential. One can then increase the 
cut-off energy until the energies for each band under consid- 
eration have converged, at which point all coefficients ci k,G ^ 
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FIG. 14: (Color online) Progressive phase update method. Starting 
from the bottom left corner of the Brillouin zone mesh (blue dots), 
the phase of the neighbor to the right (connected by the black ar- 
row) is adjusted such that Im[lnM' k,b '] = iVi/M for the pair, where 
fix i is the Berry phase in this direction. The same adjustment is 
made for the next neighbor and so on until the end of the mesh is 
reached (lower-right corner ). One then has Im[lnM( k ' b )] = v> xa /M 
for all mesh points along this path. The process is repeated for each 
path in the next reciprocal lattice direction as shown by the red ar- 
rows. 



of significant magnitude describing the Bloch periodic func- 
tions (r) are included. Typically, this requires the cut-off 
energy to be an order of magnitude greater than the upper-end 
of the energy range of interest. In our calculations a cut-off 
energy of iicut-off = 50Er is sufficient for band convergence 
and suitably limits the total number of coefficients such that 
even in three dimensions our procedure is not computationally 
expensive. Here Er is the recoil energy, defined in Eq. Q. 

The discretization of vectors in reciprocal space corre- 
sponds to considering a finite real space lattice with periodic 
boundary conditions and M primitive unit cells in each lattice 
direction. So we expect it to be valid when M is large and 
surface effects are negligible. 



Appendix B: Algorithm for reducing inter-band spread 

Our method for reducing the inter-band spread involves tak- 
ing, for each n = l,...,J — lin turn, the J — n + \ bands 
and constructing from them an «-th band that is opti- 



mally smooth in k-space, such that the most localized Wannier 
state possible may be constructed for this band. 

For each n, the algorithm, based on Ref. fl3l . proceeds as 
follows. We calculate the Hermitian matrices 



7 w 



b 



(Bl) 



where m,p run over «, . . . ,J. We then apply, for every k in 
turn, a transformation y' k ^ = l„_i <S>X^ k \ where the unitary 
diagonalises Z (k \ i.e., Z< k ) =X< k )A%W t , with A (k ) 
diagonal, reducing the spread to 



1 



The X 1 ^ are always chosen at each k such that A n is the 
largest eigenvalue of Zn*, so this spread is as small as possi- 
ble. The procedure in this paragraph is then applied repeat- 
edly until convergence is achieved. The gauge for which i2" D 
is minimized is a convergence point lfl3l . 

On occasion the above procedure can become unstable, and 



we prevent this by replacing Eq. (Bl i by an equal weighting 



of Z;„p calculated during the current and previous iteration. 
This has no effect on the locations at which the algorithm can 
converge. 



Appendix C: Algorithm for reducing intra-band spread 

For the progressive phase update method, we smooth 
the Berry connection over loops consisting of straight lines 
through the Brillouin zone, in the directions of the recipro- 
cal lattice vectors. In the reciprocal mesh representation, the 

Berry connection at each k is given by OfebIm[lnM™' b '] 
and so uniformity across a straight loop C' k b ' going through 
k' in direction b is achieved by choosing phases such that 
the projection Im[lnM™' b ^] of the connection onto this line 
is the same at each point. Specifically, since integrat- 



ing over the 



loop must 

(k,b)l 



give the Berry phase iS^k'.b) = 



L kGC (k',b) -Im[lnM^' u; ], we set Im[lnM^' b) ] = i\, (k '. b) / M at 
each point k on the loop. The loops and the order in which 
we smooth the Berry connection across them is depicted in 
Fig. M 
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